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SUMMARY 


This  report  describes  the  investigation  of  a  new  method  for  recovering 
diffraction-limited  images  through  the  turbulent  atmosphere.  The  method  con¬ 
sists  of  an  iterative  algorithm  that  constructs  an  image  from  Fourier  modulus 
data  which  is  measured  by  stellar  speckle  interferometry.  The  results  of 
this  research  indicate  that  the  method  has  the  potential  for  providing  dif¬ 
fraction-limited  images  of  earth-orbiting  satellites. 

Image  construction  experiments  were  performed  on  Fourier  modulus  data 
computer-simulated  to  include  the  effects  of  atmospheric  turbulence,  dif¬ 
fraction,  photon  (Poisson  statistics)  noise,  and  a  finite  number  of  short- 
exposure  images.  The  quality  of  the  constructed  images  was  found  to  degrade 
in  a  gradual  and  predictable  manner  as  the  signal-to-noise  ratio  decreases. 
The  rms  error  of  the  constructed  images  was  found  to  vary  approximately  as 
the  square  root  of  the  rms  error  of  the  Fourier  modulus  data.  Diffraction- 
limited  images  were  constructed  for  levels  of  photon  noise  that  would  be  ex¬ 
pected  for  imaging  satellites  through  a  1.6-meter  telescope. 

Image  construction  experiments  were  peformed  on  the  Fourier  modulus  of  a 
number  of  different  objects  of  varying  complexity.  Interpretation  of  the 
results  was  complicated  by  a  tendency  of  the  algorithm  to  stagnate  at  local 
minimum  having  the  appearance  of  a  good  quality  image  superimposed  by  a  pat¬ 
tern  of  stripes.  Nevertheless,  the  results  were  suggestive  that  the  solution 
is  usually  unique;  only  for  an  object  satisfying  special  conditions  is  the 
Fourier  modulus  ambiguous.  As  the  signal-to-noise  ratio  decreases,  the  am¬ 
biguity  of  the  solution  increases,  but  that  ambiguity  takes  the  form  of  noise 
in  the  constructed  image  rather  than  a  complete  change  in  the  basic  shape  of 
the  image. 

Several  variations  of  the  algorithm  were  compared,  and  one  was  found  to 
converge  considerably  faster  than  the  others.  The  convergence  time  was 
longer  for  objects  of  greater  complexity.  For  an  array  size  of  128  x  128 
pixels,  typically  one  hundred  iterations  are  required,  taking  about  two 
minutes  on  an  array  processor. 
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FOURIER  MODULUS  IMAGE  CONSTRUCTION 


1 

INTRODUCTION 

The  purpose  of  this  research  project  is  to  investigate  the  effectiveness 
of  a  new  method  of  recovering  fine  resolution  imagery  through  the  turbulent 
atmosphere. 

Ordinarily,  atmospheric  turbulence  limits  the  resolution  of  an  image  ob¬ 
tained  through  a  large  optical  telescope  to  about  one  second  of  arc  at  best, 
which  corresponds  to  the  diffraction  limit  of  an  optical  telescope  of  aper¬ 
ture  10  cm.  By  comparison,  the  theoretical  diffraction-limited  resolution 
of  a  2-meter  telescope  with  no  atmospheric  turbulence  is  0.05  arc-seconds  — 
20  times  finer.  Compensated  imaging  (real-time  correction  of  atmospheric 
turbulence-induced  phase  errors  by  adaptive  optics)  can  produce  diffraction- 
limited  images,  but  the  system  is  expensive.  An  inexpensive  alternative  ap¬ 
proach  is  to  gather  fine  resolution  information  by  stellar  speckle  inter¬ 
ferometry  [1],  and  then  process  that  information  to  form  an  image.  The  in¬ 
formation  obtained  from  speckle  interferometry  is  the  modulus  of  the  object's 
Fourier  transform,  or  equivalently  the  autocorrelation  of  the  object.  An 
image  can  be  constructed  from  this  information  using  an  iterative  algorithm 
[2].  The  imagery  thus  obtained  would  have  diffraction-limited  resolution. 

This  research  project  had  three  major  goals:  (1)  to  determine  how  the 
quality  of  the  constructed  images  varied  with  the  signal-to-noise  ratio  of 
the  Fourier  modulus  data;  (2)  to  determine  whether  the  algorithm  could  pro¬ 
duce  spurious  results  (that  is,  whether  the  constructed  image  is  unique); 
and  (3)  to  investigate  changes  in  the  algorithm  to  improve  convergence  of 
the  algorithm  and  to  determine  how  the  convergence  time  varies  with  the  pro¬ 
perties  of  the  Fourier  modulus  data.  Section  2  of  this  report  gives  back¬ 
ground  to  the  iterative  algorithm.  Sections  3  through  5,  respectively,  de¬ 
scribe  the  results  relating  to  the  three  goals  listed  above.  Section  6  sum¬ 
marizes  the  conclusions  and  outlines  areas  where  further  work  is  needed. 
Appendix  A  contains  a  proof  of  convergence  for  one  version  of  the  iterative 
al gor ithm. 
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2 

ITERATIVE  METHOD  BACKGROUND 


This  section  briefly  describes  the  iterative  method,  which  is  discussed 
in  more  detail  elsewhere  [2,  3]  and  in  Section  5.  Since  the  experiments  were 
performed  on  a  digital  computer,  the  object  and  its  Fourier  transform  are 
sampled  functions.  Throughout  this  report,  the  convention  is  used  that 
functions  represented  by  upper  case  letters  are  the  Fourier  transforms  of 
functions  represented  by  lower  case  letters. 

Let  the  object  distribution  be  f  and  its  discrete  Fourier  transform 
be 


Fpq  s  |Fpq|  exP(iePq)  “^fmn]  "  2  fmn  exp^ (2*/N) (mp  ♦  nq)]  (1) 

mn 


where  m,  n,  p,  and  q  »  0,  1 ,  2,  .  .  N  -  1 .  m  and  n  are  referred  to  as 

spatial  coordinates  and  p  and  q  as  spatial  frequency  coordinates.  The  object 

brightness  distribution  f  is  real  and  nonnegative  and  F...  Is  Hermi- 

mn  pq 

tian.  We  assume  that  only  |F  |  is  measured,  although  in  some  cases  of 

practical  interest  «pq  may  be  measurable  for  very  low  spatial  frequen¬ 

cies.  The  diameter  of  the  object  can  be  computed  since  it  is  half  the  diam¬ 
eter  of  the  autocorrelation  f  ★  f  1  [ |  Fpq |^3.  The  problem  is  to 

reconstruct  the  object  fmn,  or  equivalently  the  phase  epq,  from  j  Fpq  | 

consistent  with  the  fact  that  f  t  0.  Use  of  the  constraint  of  the  ob- 

mn 

ject's  computed  diameter  is  not  essential,  since  any  nonnegative  solution 
having  a  Fourier  modulus  j  Fqq  j  automatically  has  the  correct  diameter. 

A  block  diagram  of  the  iterative  method  is  depicted  in  Figure  la.  The 
four  steps  of  each  iteration  are  as  follows.  (1)  An  input  image  gmn 
(which  in  some  instances  may  be  considered  an  estimate  of  the  object)  is 
Fourier  transformed: 

G__  =  I  G_J  exp( id  )  %<jT[gm„] 
pq  |  pq|  rv  pq'  nt^L:’mnJ 
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Figure  1.  (A)  Block  diagram  of  the  Iterative  method; 

(B)  the  Input-output  concept  for  the  iterative  method. 
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(2)  The  result  is  made  to  satisfy  the  measured  Fourier  modulus  and  the  com¬ 
puted  phase  is  left  unchanged: 


exP(i*pq) 


(3)  The  result  is  inverse  Fourier  transform: 


Vm 

(4)  Based  on  how  g'mn  violates  the  object-domain  constraints,  a  new  input 
gmn  is  formed,  and  the  process  is  repeated. 

Taken  together  as  a  group,  the  first  three  steps  can  be  viewed  as  a  non¬ 
linear  system  having  an  Input  gmn  and  an  output  g'mn  as  depicted  in 
Figure  lb.  The  output  g'mn  is  guaranteed  to  have  a  Fourier  transform  with 
the  measured  modulus.  Therefore,  if  g'mn  is  nonnegative,  then  it  is  a  so¬ 
lution  to  the  problem;  that  is,  it  simultaneously  agrees  with  the  measured 
data  in  the  Fourier  domain  and  the  a  priori  nonnegativity  constraint  in  the 
image  (or  object)  domain.  The  procedure  is  to  iteratively  change  the  input 
in  such  a  way  as  to  drive  the  output  to  be  everywhere  nonnegative.  Driving 
the  output  to  be  zero  outside  the  computed  diameter  may  also  be  Included; 
imposition  of  the  diameter  constraint  is  not  essential,  but  it  decreases  the 
convergence  time  of  the  algorithm. 

The  iterative  method  embraces  a  family  of  algorithms  for  altering  the 
input  in  order  to  drive  the  output  to  satisfy  the  object-domain  constraints. 
Two  such  methods  of  altering  the  input  that  were  used  for  the  experiments 
reported  here  are  given  by 


and 


9k+l ,mn 


3k,mn-  * 

(  0  ,  (m,n)  e  y 


Vi, 


mn 


9k, m  ■  (”•")*  * 


9k, mn  ‘  89k,mn-  <m-">  *  * 


(2) 


(3) 


12 


where  g.,  and  g*.  are  the  input  and  output,  respectively,  at  the 
ktn  iteration,  y  is  the  set  of  points  at  which  g'^  mn  violates  the  con¬ 
straints,  and  0  is  a  constant.  We  refer  to  the  iterative  method  using  Eq. 
(2)  as  the  error-reduction  approach,  since  it  can  be  shown  that  the  mean- 
squared  error  defined  in  the  Fourier  domain  by 


or  in  the  image  domain  by 


can  only  decrease  at  each  iteration.  The  proof  of  convergence  is  given  in 
Appendix  A.  These  two  error  measures  indicate  the  degree  to  which  the  solu¬ 
tion  agrees  with  the  measured  Fourier  modulus.  Unfortunately,  the  fractional 
decrease  of  the  error  becomes  very  small  as  the  iterations  progress,  and  the 
error-reduction  approach  usually  does  not  converge  to  a  solution  in  a  prac¬ 
tical  sense. 

An  important  notion  that  is  used  in  Eq.  (3)  and  other  related  methods  of 
altering  the  input  is  that  a  change  of  the  input  tends  to  result  in  a  similar 
(but  somewhat  different)  change  of  the  output  [41.  We  refer  to  these  methods 
of  altering  the  input  as  the  input-output  approach.  Although  the  mean- 
squared  error  does  not  necessarily  decrease  at  each  iteration  (it  may  even 
increase),  in  practice  the  input-output  approach  converges  much  faster  than 
the  error-reduction  approach.  The  optimum  choice  of  the  constant  s  in  Eq. 
(3)  depends  on  the  statistics  of  and  F„„.  We  have  found  values  of  e 

pq  pq 

between  0.5  and  1.0  to  work  well.  The  input-output  approach  can  be  used  by 
itself,  or  it  can  be  alternated  with  the  error-reduction  approach.  For  the 
results  shown  in  Sections  3  and  4,  we  alternated  between  Eq.  (2)  and  Eq.  (3), 
typically  performing  thirty  to  sixty  iterations  using  Eq.  (3)  followed  by 
ten  iterations  using  Eq.  (2). 
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The  region  y  in  Eqs.  (2)  and  (3)  includes  all  points  at  which  g’k 
violates  the  object-domain  constraints.  The  primary  constraint  is  that  the 
image  be ’ nonnegative.  Another  constraint  that  may  be  imposed  is  that  the 
diameter  of  the  image  not  exceed  half  the  diameter  of  the  autocorrelation. 
For  the  set  of  experiments  described  in  Sections  3  and  4,  we  applied  a  tight 
diameter  constraint  for  the  first  few  iterations,  but  only  constrained  the 
image  to  half  the  field-of-vlew  (and  the  Images  always  fit  within  less  than 
half  the  field-of-vlew)  in  each  dimension  for  the  bulk  of  the  Iterations. 


As  mentioned  earlier,  the  diameter  constraint  is  not  essential,  since 
all  real,  nonnegative  solutions  automatically  have  the  correct  diameter  (half 
the  diameter  of  the  autocorrelation).  However,  during  the  iterations  the 
intermediate  results  (the  output  images  g')  do  have  a  diameter  exceeding  the 
known  diameter,  and  by  applying  the  diameter  constraint  the  convergence  time 
of  the  algorithm  is  decreased.  That  g'  can  have  a  diameter  greater  than  half 
the  diameter  of  the  autocorrelation  can  be  seen  from  the  following  example. 
Suppose  that  for  a  given  iteration  the  input  g  is  nonnegative  and  has  the 
correct  diameter,  but  has  an  incorrect  Fourier  modulus.  After  Fourier  trans¬ 
formation,  its  Fourier  transform  G  is  modified  to  satisfy  the  measured 
Fourier  modulus,  |F  |,  which  we  will  assume  is  exactly  the  Fourier  modulus  of 
the  object.  Neglecting  the  case  where  |  G|  *  0  (it  is  almost  never  exactly 
equal  to  zero,  due  to  noise,  etc),  this  can  be  accomplished  by  multiplying  G 

by  I F | / | G | , 


I 


pq! 


By  the  convolution  theorem  the  resulting  output  image  can  be  expressed  as 


^mn 


In  general  the  Fourier  transform  of  (|F|/|G|)  will  have  positive  and  negative 
sidelobes  that  extend  over  the  entire  field-of-view.  Therefore  g'  will 
contain  areas  of  positive  and  negative  values  that  extend  beyond  half  the 
diameter  of  the  autocorrelation.  Nevertheless  g1  and  the  object  must  have 
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identical  autocorrelations:  since  |G'|  «  | F|  ,  g' 

^[|F|^]  ■  f  ★  f.  The  diameter  of  g'  can  exceed  half  the  diameter  of 
the  autocorrelation  only  if  g'  has  regions  of  negative  value,  since  negative 
values  are  required  to  get  the  cancellation  needed  to  reduce  the  diameter  of 
the  autocorrelation  to  less  than  twice  the  diameter  of  g'. 


3 

NOISE  PROPERTIES  OF  IMAGES  CONSTRUCTED  FROM  FOURIER  MODULUS 


This  section  describes  a  set  of  experiments  aimed  at  determining  how  the 
quality  of  the  constructed  images  depends  on  the  amount  of  noise  present  in 
the  'ured  Fourier  modulus  data.  Noisy  Fourier  modulus  data  was  simulated 
having  predetermined  amounts  of  noise,  images  were  constructed  from  the  noisy 
Fourier  modulus  data,  and  a  measure  was  made  of  the  quality  of  the  con¬ 
structed  images. 


3.1  SIMULATION  OF  NOISY  FOURIER  MODULUS  DATA 


The  Fourier  modulus  data  used  for  the  reconstruction  experiments  was 
simulated  to  include  the  effects  of  atmospheric  turbulence  and  photon  noise. 
The  object  used  (a  digitized  photograph  of  a  P72-2  sensor  testbed)  is  shown 
in  Figure  2,  as  it  would  appear  through  a  diffraction-limited  telescope  with 
no  atmospheric  turbulence  and  no  noise.  The  object  is  about  64  x  40  pixels 
in  extent,  imbedded  in  a  128  x  128  array.  The  object  f  was  convolved 
with  156  different  point-spread  functions  (PSFs)  to  produce  156  different 
blurred  images 


th 


f  *  i>nn  sisn. 
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156  (6) 


is  the  kL"  PSF.  The  object  and  the  PSFs  were  supplied  by 

mn 

B.  L.  McGlamery  (Visibility  Laboratory,  Scripps  Institution  of  Oceanography, 
U.C.  San  Diego)  and  were  computed  from  phase  functions  representing  the  ef¬ 
fects  of  a  model  of  atmospheric  turbulence  having  a  Kolmogorov  spectrum  [5], 

assuming  a  telescope  lens  diameter  of  1.6  meters  and  rn  (Fried's  parameter 

0  Ik  1 

[6])  =  0.1  meters.  Examples  of  the  blurred  images  d^n  are  shown  in 

Figure  3(a)— (d).  The  blurred  images  were  then  converted  to  mean  photon 
counts  and  subjected  to  a  Poisson  noise  process:  each  pixel  was  replaced 

with  a  sample  drawn  from  a  Poisson  distribution  with  mean  and  variance  equal 


to  the  mean  photon  count,  resulting  in  the  noisy  degraded  images  i 


(k) 

mn 


A  point  of  reference  for  the  amount  of  photon  noise  present  in  the  data 
is  given  by  an  analysis  done  by  B.  L.  McGlamery  [7],  He  assumed  that  the 
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Fiyure  2.  The  reference  object. 


!  ?  t% 


Figure  3.  Degraded  images.  (A)-(D)  Noise-free  blurred  images;  degraded 
(blurred  and  noisy)  images —number  of  photons  per  degraded  images: 
(E)-(H)  305,000;  (I)-(L)  6143;  (M)-(P)  643. 


short-exposure  blurred  images  were  obtained  with  a  system  having  the  follow¬ 
ing  parameters: 


1.6  m 

=  telescope  diameter 

5  msec 

■  exposure  time 

100  nm 

-  spectral  bandwidth 

550  nm 

-  center  wavelength 

0.8 

■  atmospheric  transmittance 

and  assuming  a  5900*K  spectral  curve  for  the  solar  Irradlance,  and  an  S20 
sensor  spectral  sensitivity  curve.  For  the  P72-2  object  (Figure  2)  at  a  mag¬ 
nification  equivalent  to  a  distance  of  300  km,  the  blurred  Images  (reflec¬ 
tance  maps)  were  multiplied  by  a  constant  of  value  2779  to  convert  the  Image 
to  a  mean  value  of  photoelectrons/pixels  (mean  photon  counts).  We  refer  to 
this  constant  for  converting  reflectance  Into  mean  photoelectrons  as  the 
photoelectron  scale  factor.  The  value  of  2779  pertains  to  the  system  and 
exposure  time  described  above,  which  was  considered  to  be  realistic  for  ob¬ 
serving  satellites  (although  the  spectral  bandwidth  of  100  nm  may  be  con¬ 
sidered  to  be  too  high  by  as  much  as  a  factor  of  two).  Using  156  blurred 
Images,  this  case  serves  as  a  baseline.  For  other  cases  there  would  be 
greater  or  fewer  photons  available  according  to  the  reflectivity  of  the  ob¬ 
ject  and  the  values  of  the  parameters  listed  above.  The  net  effect  on  the 
photon  noise  of  all  of  these  factors  was  simulated  by  lumping  them  all  into 
one  number,  the  photoelectron  scale  factor.  For  our  experiments  the  photo¬ 
electron  scale  factor  was  varied  between  6  and  55580  which  is  equivalent  to 
between  0.002  times  and  20  times  the  number  of  photons  for  the  baseline  case. 
The  values  of  the  photoelectron  scale  factor  that  were  used  are  55580,  27790, 
13895,  5558,  2779,  1390,  556,  278,  139,  56,  28,  14,  and  6. 

Examples  of  the  noisy,  degraded  images,  including  both  the  blurring  due 
to  the  atmosphere  (and  the  telescope  aperture)  and  photon  noise,  are  shown 
in  Figure  3. 

This  simulated  speckle  data  was  then  processed  by  the  Labeyrie  method 
[1],  as  modified  by  Goodman  and  Belsher  [8]  and  later  by  others  [9],  to  ar¬ 
rive  at  an  estimate  of  the  modulus  of  the  Fourier  transform  of  the  object: 
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where 

Nc  -EE  W  -p'S’  (8) 

is  a  constant  equal  to  the  total  number  of  photons  detected  [8,  9], 

E,S^|2  is  the  square  of  the  MTF  of  the  speckle  interferometry  pro¬ 
cess,  and  1^  is  a  weighting  factor.  Examples  of  |Fp(J  are  shown  in 
Figure  4.  The  subtraction  of  Nc  eliminates  one  particular  constant  noise 
term  due  to  photon  statistics.  Other  noise  terms  remain.  If  the  sensor  is 
operating  in  a  photon-counting  mode,  then  Nc  can  be  calculated  from  Eq. 

(8)  as  it  was  for  our  experiments;  otherwise,  Nc  can  be  estimated  as  the 
minimum  value  of  \  at  hi9h  spatial  frequencies.  In  a  real- 

world  situation,  it  may  also  be  necessary  to  compensate  for  the  MTF  due  to 
the  detection  process  before  subtracting  Nc  [10].  In  the  real  world  an 

estimate  of  the  MTF  in  the  denominator  of  Eq.  (7)  can  be  obtained  from  mea¬ 

surements  on  an  unresolved  star  or  determined  from  a  model  of  the  atmosphere. 
For  this  experiment  the  intent  was  to  separate  out  the  effects  due  to  photon 
noise  alone,  so  the  exact  MTF  was  used.  Division  by  the  MTF  also  compensates 
for  the  MTF  due  to  the  telescope  aperture,  resulting  in  a  net  system  transfer 
function  that  is  unity  within  a  circle  of  radius  equal  to  twice  the  radius 
of  the  telescope  aperture,  and  is  set  to  zero  outside  that  aperture.  The 
result  would  be  a  system  impulse  response  having  undesirable  negative  side- 
lobes.  The  weighting  function  WpC|  was  included  in  Eq.  (7)  in  order  to  re¬ 
store  the  natural  MTF  associated  with  a  telescope  having  a  circular  aperture 
(that  MTF  being  the  autocorrelation  of  the  circular  aperture),  and  the  re¬ 
sulting  system  impulse  response  no  longer  has  negative  sidelobes.  The  unde¬ 
graded  image  shown  in  Figure  2,  which  we  refer  to  as  the  reference  object, 
was  also  subjected  to  the  same  weighting  function  in  the  Fourier  domain  to 
impose  on  it  the  transfer  function  of  the  telescope  aperture.  The  radius  of 
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Figure  4.  Fourier  modulus  estimates  computed  from  156  degraded  images  each. 
Number  of  photons  per  degraded  image:  (A)  infinity  (noi se-free) ;  (B)  305,000; 

(C)  6143;  (D)  643. 


Wpq  was  62  pixels,  corresponding  to  the  MTF  of  a  telescope  of  diameter  1.6 
meters. 

An  indication  of  the  signal-to-noise  ratio  of  the  Fourier  modulus  esti- 
A 

mate  |F  |  can  be  obtained  from  the  total  number  of  photons  detected, 

N  .  Values  of  N  for  the  Fourier  modulus  estimates  were  varied  from  9.5 
c  8  c 

x  10  for  the  lowest  noise  case  (photoelectron  scale  factor  =  55580)  to 

5 

i.O  x  10  for  the  highest  noise  case  (photoelectron  scale  factor  =  6). 
Since  156  degraded  images  were  used  to  compute  the  Fourier  modulus  esti¬ 
mates,  this  corresponds  to  a  range  of  6.1  x  10®  to  643  photons  per  de¬ 
graded  image.  Since  each  blurred  image  extended  over  an  area  of  about  64  by 
64  pixels,  this  corresponds  to  a  range  of  1490  to  0.16  photons  per  pixel  per 
degraded  image  over  the  area  of  the  degraded  image.  For  the  baseline  case 

having  a  photoelectron  scale  factor  of  2779,  the  value  of  N  was  4.7  x 

7  5  c 

10  ,  which  is  equivalent  to  3  x  10  photons  per  degraded  image,  or  74 

photons  per  pixel  per  degraded  image  over  the  area  of  the  degraded  image. 

Although  it  would  not  be  measurable  from  telescope  data,  for  our  computer 

simulations  we  also  computed  the  actual  mean-squared  error  of  the  Fourier 

modulus  estimate: 


where  Fp^  is  assumed  to  include  the  weighting  factor  Wp^.  A  plot  of 
Ej^j,  the  rms  error  of  the  Fourier  modulus  estimate,  versus  the  number  of 
photons  per  degraded  image  is  shown  in  the  upper  curve  of  Figure  5. 

Since  the  signal-to-noise  ratio  of  the  Fourier  modulus  estimate  tends  to 
decrease  with  increasing  spatial  frequencies,  one  would  hope  to  reconstruct 
low-frequency  (low  resolution)  features  of  an  object  even  under  very  noisy 
conditions  for  which  the  high-frequency  information  is  lost.  For  this  rea¬ 
son,  reconstruction  experiments  were  also  performed  on  low  resolution  ver¬ 
sions  of  the  Fourier  modulus  estimates.  Formation  of  low  resolution  versions 
of  |^|  were  accomplished  by  replacing  the  radius-62  in  Eq.  (7)  by  a 
Wpg  of  the  same  form,  but  of  radius  16  pixels  ( i . e . ,  approximately  four 
times  coarser  resolution  than  the  baseline  "high  resolution"  version).  The 
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NO  OF  PHOTONS  PER  DEGRADED  IMAGE  (Nc/156) 


Figure  5.  Rms  error  of  Fourier  modulus  estimate,  EijTi,  versus  number  of 
photons  per  degraded  Image,  Nc/156.  Upper  curve:  high  resolution  version; 
lower  curve:  low  resolution  version. 
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lower  curve  in  Figure  5  shows  EiA  versus  the  number  of  photons  per  de- 

!AI 

graded  mage  for  the  low  resolution  version  of  |F  |.  The  fact  that  the  curve 
for  the  low-resolution  case  is  below  that  of  the  high-resolution  case  (that 
is,  the  error  of  the  Fourier  modulus  estimate  is  lower)  confirms  the  notion 
that  the  error  increases  with  increasing  spatial  frequencies. 


3.2  IMAGE  CONSTRUCTION  RESULTS  FROM  NOISY  FOURIER  MODULUS 

Image  construction  by  the  iterative  method  (as  described  in  Section  21 
was  carried  out  on  13  high-resolution  and  8  low-resolution  Fourier  modulus 
estimates  which  were  generated  in  the  manner  described  in  Section  3.1.  In 
each  case,  although  useful  images  were  usually  available  after  about  100 
iterations,  600  iterations  (the  same  sequence  of  iterations  for  all  cases ) 
were  carried  out  in  order  to  ensure  that  the  convergence  of  the  algorithm 
was  as  complete  as  was  practical. 

The  rms  error  E  .  which  is  a  measure  of  how  well  the  constructed  imaqe 

agrees  with  the  object-domain  constraints  and  Fourier  modulus  estimate,  was 

reduced  to  about  0.007  for  the  lowest-noise  case,  to  0.05  for  the  median- 

noise  case,  and  to  0.4  for  the  highest-noise  case.  In  no  case  was  EQ 

driven  to  zero;  that  is,  in  no  case  did  the  algorithm  converge  to  a  solution 

that  was  in  perfect  agreement  both  with  the  Fourier  modulus  estimate  and  with 

the  object’s  nonnegativity  constraint.  Such  agreement  would  be  impossible 

because  the  noise  present  in  the  Fourier  modulus  estimate  makes  that  estimate 

inconsistent  with  the  object's  nonnegativity  constraint.  In  particular,  if 

the  corresponding  estimate  of  the  object's  autocorrelation  is  computed  by 

,  A ,? 

Fourier  transforming  |  F  |  ,  that  autocorrelation  estimate  is  found  to  have 
regions  of  negative  values.  It  can  easily  be  shown  that  an  autocorrelation 
with  negative  values  can  arise  only  from  an  object  with  negative  values. 
That  is,  there  can  be  no  nonnegative  (physical)  object  that  could  give  rise 
to  the  Fourier  modulus  estimate  1^1.  (In  addition,  even  if  the  estimate  of 
the  object's  autocorrelation  were  nonnegative,  then  it  is  still  possible  that 
there  would  be  no  nonnegative  object  giving  rise  to  it.  Although  we  do  not 
know  in  general  how  to  tell  if  a  given  function  is  an  autocorrelation  func¬ 
tion,  it  can  be  shown  that  not  all  nonnegative  symmetric  functions  are  auto- 
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correlation  functions,  and  a  sum  of  autocorrelation  functions  is  not  neces¬ 
sarily  an  autocorrelation  function.)  Nevertheless,  this  iterative  image 
construction  method,  which  relies  on  the  nonnegativity  of  the  solution  and 
strives  toward  a  nonnegative  solution,  finds  a  solution  that  has  a  minimum 
amount  of  negative  values;  in  doing  so,  it  constructs  an  image  that  is  (at 
least  for  cases  where  only  a  moderate  amount  of  noise  is  present)  a  good  ap¬ 
proximation  to  the  original  object. 

The  image  construction  results  are  shown  in  Figure  6  and  7  for  the  high- 
resolution  and  low-resolution  cases,  respectively.  Images  that  came  out  in¬ 
verted  (rotated  in  the  plane  of  the  page  by  180*)  were  re-inverted  in  order 
to  take  on  the  same  orientation  as  the  reference  object.  Inverted  solutions 
are  always  possible  since  the  Fourier  modulus  of  f ( — x )  equals  the  Fourier 
modulus  of  f(x),  for  f(x)  real  valued.  Figure  6(E)  shows  the  result  for  the 
baseline  case.  In  the  high-resolution  case,  a  good  quality  image,  Figure 
6(H),  was  constructed  for  a  photon  count  (slgnal-to-noise  ratio)  ten  times 
worse  than  the  baseline  case.  A  degraded  but  still  recognizable  image,  Fig¬ 
ure  6(K),  was  constructed  for  a  photon  count  one  hundred  times  worse  than 
the  baseline  case.  At  four  hundred  times  worse  photon  count  than  the  base¬ 
line  case,  Figure  6(H),  it  appears  that  no  useful  image  information  remains. 

Also  shown  for  comparison  in  Figures  6(N)-6(U)  are  short-exposure  images 
with  no  atmosphere  present  (diffraction-limited  and  photon- 1 imi ted)  for  the 
eight  lowest  signal-to-noise  ratio  cases,  corresponding  to  Figures  6 ( F )— 6 (M ) , 
respectively.  From  this  it  is  seen  that  a  single  short-exposure  image  with 
no  atmosphere  present  is  better  than  an  image  constructed  from  the  Fourier 
modulus  estimate  based  on  156  short-exposure  images.  One  might  also  wish  to 
compare  the  image  construction  results  with  the  sum  of  156  photon-limited 
images  with  no  atmosphere  present.  For  the  lowest  signal-to-noise  ratio  case 
having  photoelectron  scale  factor  of  6,  the  sum  of  156  images  would  be  the 
same  as  a  single  image  of  photelectron  scale  factor  of  936  =  6  x  156,  which 
corresponds  to  a  result  half  way  between  Figures  6(N)  and  6(0). 

The  low  resolution  results  shown  in  Figure  7  are  for  the  eight  lowest 
signal-to-noise  ratio  cases  only.  Comparing  Figure  7  with  Figure  6(F)-6(M), 
for  the  same  number  of  photon  counts,  N  ,  the  constructed  images  appear  to 
be  less  noisy  in  the  low-resolution  case  than  in  the  high-resolution  case. 
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Figure  6.  (A)  Reference  object;  (B)-(M)  images  constructed  from  Fourier 

modulus  estimates,  corresponding  to  the  data  points  in  Figure  5  for  the 
high  resolution  case,  excluding  the  lowest  noise  case.  The  baseline  case 
is  (E).  Each  third  image  corresponds  to  a  factor  of  10  in  number  of 

photons  detected. 
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Figure  7.  (A)  Low  resolution  version  of  reference  object;  (B)-(I)  images 

constructed  from  low  resolution  versions  of  the  Fourier  modulus  estimates, 
corresponding  to  the  data  points  in  the  lower  curve  of  Figure  5. 
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It  should  be  noted,  however,  that  in  the  high-resolution  case,  the  noise  in 
the  constructed  image  is  mostly  in  the  higher-frequency  details,  and  the 
lower-frequency  components  (the  coarser  details)  of  the  constructed  images 
remain  more  faithful  than  the  higher-frequency  details.  Based  on  just  their 
low  frequency  content  (coarse  details),  the  high-resolution  and  low- 
resolution  cases  appear  to  be  comparable. 

Although  It  would  not  be  measurable  from  telescope  data,  for  our  digital 
simulations  we  computed  the  actual  mean-squared  error  of  the  constructed  Im¬ 
age,  since  we  have  available  the  reference  object.  Before  this  mean-squared 
error  can  be  computed,  the  constructed  Image  must  be  brought  into  coincidence 
with  the  reference  object,  since  the  constructed  Image  may  be  translated  and 
even  spatially  Inverted  relative  to  the  reference  object  [2],  The  choice  of 
whether  to  invert  the  image  and  what  translation  Is  required  was  determined 
by  successively  computing  the  cross-correlations  of  the  reference  object  with 
the  constructed  image  at  normal  and  inverted  orientations,  and  searching  for 
the  maximum  over  the  two  cross-correlations.  Inversion  of  the  image  (when 
necessary)  was  effected  by  conjugating  its  Fourier  transform,  and  translation 
of  the  image  was  accomplished  by  multiplying  its  Fourier  transform  by  a  lin¬ 
ear  phase  factor.  The  resulting  constructed  complex  Fourier  transform 
was  used  to  compute  the  mean-squared  error 
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which  is  equivalent  to  the  mean-squared  error  defined  in  the  image  domain  as 
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where  gmr)  is  the  translated  and  possibly  inverted  version  of  the  con¬ 
structed  image.  In  Eqs.  (10)  and  (11)  the  reference  object  f  and  its 
Fourier  transform  are  assumed  to  have  been  subjected  to  the  same  weighting 
function  Vi  as  the  Fourier  modulus  estimate. 
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The  rms  error  of  the  constructed  images,  is  shown  as  a  function  of 
E|£|,  the  rms  error  of  the  Fourier  modulus  estimate,  in  Figure  8,  which 
shows  that  E  appears  to  vary  approximately  as  the  square  root  of  E|£|. 
This  approximate  square  root  relationship  was  not  expected  —  no  theoretical 
results  are  known  that  predict  this  relationship.  For  a  given  Fj^j,  E  is 
approximately  the  same  for  the  low  resolution  cases  (triangles)  as  for  the 
high  resolution  cases  (circles).  Therefore,  except  for  possible  savings  in 
computation  time  there  does  not  seem  to  be  any  advantage  to  low-pass  filter¬ 
ing  the  data  before  applying  the  iterative  algorithm.  A  more  sophisticated 
filter  such  as  Wiener  filter  could,  on  the  other  hand,  prove  advantageous. 

For  the  low  noise  cases  E  did  not  go  below  0.1  although  the  square  root 
of  E|^|  was  well  below  that  level.  The  cause  of  this  error  in  the  low 
noise  cases  seems  to  have  resulted  not  from  the  noise  but  from  the  stagnation 
of  the  algorithm  at  an  image  that  perhaps  is  at  a  local  minimum  of  error, 
but  not  at  a  global  minimum.  These  images  look  like  the  reference  object 
with  a  low  contrast  sinusoidal  pattern  superimposed  (or  as  stripes  across 
the  image).  Figures  9 ( a )— 9 (d )  show  four  different  images  constructed  from 
the  same  Fourier  modulus  estimate  (Nc  =  9.5  x  10  ),  each  arrived  at  by 
the  same  sequence  of  iterations  but  with  different  arrays  of  random  numbers 
as  the  initial  input  to  the  first  iteration.  The  stripes  appear  in  each  im¬ 
age,  but  with  different  orientations  and  spatial  frequencies.  This  pheno¬ 
menon  was  noticed  earlier  [2],  but  a  complete  explanation  of  it  has  not  been 
available.  Further  discussion  of  this  phenomenon  is  given  in  Section  5. 
Also  shown  in  Figure  9  are  four  different  images  constructed  from  the  same 
Fourier  modulus  for  a  higher  noise  case  (Nc  =  9.5  x  10®).  In  this  case 
the  stripes  do  not  occur;  but  a  noticeable  amount  of  random  noise  occurs  in 
the  background  of  the  images,  and  the  details  of  that  noise  are  different 
for  each  constructed  image. 

The  image  construction  results  shown  in  Figure  6  were  arrived  at  from 
initial  inputs  consisting  of  random  numbers.  In  addition,  images  were  also 
constructed  usinq  the  reference  object  as  the  initial  input.  Although  such 
an  initial  input  is  unrealistic  from  the  point  of  view  that  in  the  real  world 
the  reference  object  would  not  be  known,  those  image  construction  results 
were  interesting  because  they  shed  some  light  on  the  problem  of  the  stripes. 
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RMS  FOURIER  MODULUS  ERROR 

Figure  8.  Rms  error  of  constructed  Images,  E,  versus  rms  error  of  the 
Fourier  modulus  estimates,  Ejpj .  Circles  Indicate  high  resolution  versions 

and  triangles  Indicate  low  resolution  versions.  Squares  Indicate  results 
when  the  reference  object  was  used  as  the  Initial  Input  to  the  algorithm 

(see  text).  The  solid  line  shows  Y  pi  comparison. 
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Figure  9.  Four  images  constructed  using  different  starting  inputs. 
(A)-(D)  Low  noise  case  (6.1  x  106  photons  per  degraded  image)  showing 
stripes;  (E)-(H)  higher  noise  case  (6.1  x  l(r  per  photons  per  degraded 
image).  These  pictures  were  intentionally  overexposed  in  order  to  em¬ 
phasize  the  stripes  and  the  background  noise. 


When  starting  with  the  reference  object  as  the  initial  input,  we  have  some¬ 
thing  that  is  not  in  perfect  agreement  with  the  noisy  Fourier  modulus  esti¬ 
mate.  As  the  iterations  progress,  the  output  image  departs  from  tl>e  refer¬ 
ence  object  while  becoming  in  better  agreement  with  the  Fourier  modulus  esti¬ 
mate.  Finally,  the  iterations  stagnate  at  a  point  where  the  result  is  in  as 
good  agreement  as  possible  with  the  simultaneous  constraints  of  object  non¬ 
negativity  and  Fourier  modulus  equalling  the  Fourier  modulus  estimate. 

For  the  medium  to  high  noise  cases,  the  constructed  Images  resulting  from 
the  reference  object  as  the  Initial  Input  were  comparable  to  the  constructed 
images  resulting  from  random  numbers  as  the  Initial  Input;  In  addition,  the 
error  measures  EQ  and  ?  were  comparable  for  the  two  cases.  For  the  low 
noise  cases,  the  constructed  Images  resulting  from  the  reference  object  as 
the  initial  Input  were  comparable  In  appearance  to  the  constructed  Images 
resulting  from  random  numbers  as  the  Initial  Input  except  that  no  stripes 
were  present.  In  addition,  both  the  error  measures  EQ  and  E  were  consider¬ 
ably  less  than  the  values  for  the  case  of  the  random  Initial  Input.  E  for 

this  case  Is  plotted  In  Figure  8  along  with  the  case  of  the  random  initial 

Input,  where  It  Is  seen  that  t  approaches  zero  as  E|£|  approaches  zero,  as 
one  would  expect.  For  the  four  lowest  noise  cases,  the  error  EQ  for  the 

case  of  the  random  Initial  input  versus  the  case  of  the  reference  object  as 
initial  input  was  0.00737  vs  0.00333,  0.00748  vs  0.00391,  0.00794  vs  0.00526, 
and  0.0114  vs  0.00933.  For  higher  amounts  of  noise  EQ  was  comparable  in 
the  two  cases.  Therefore,  for  the  low  noise  cases  one  can  distinguish  the 
local  minimum  having  stripes  from  the  global  minimum  not  having  stripes  by 

the  value  of  the  rms  error  E„.  It  Is  this  fact  that  a  smaller  error  E. 

o  o 

occurred  for  the  constructed  Image  not  having  stripes  that  leads  us  to  be¬ 
lieve  that  the  stripes  represent  a  problem  of  encountering  a  local  minimum 
and  do  not  represent  a  fundamental  uniqueness  problem. 

The  most  Important  result  from  this  set  of  experiments  on  the  noise  pro¬ 
perties  is  that  as  the  noise  increases,  the  quality  of  the  constructed  images 
degrades  in  a  gradual  and  predictable  manner.  Furthermore,  the  quality  of 
the  constructed  images  Is  very  good  for  the  realistic  slgnal-to-noise  ratios 
expected  for  Imaging  satellites.  Images  of  this  quality  were  obtained  in 
spite  of  the  fact  that  Wiener  filtering  was  not  performed.  By  adding  a  post 
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processing  step  of  Wiener  filtering  the  images,  improved  image  quality  may 
be  expected. 
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4 

UNIQUENESS  OF  IMAGES  CONSTRUCTED  FROM  FOURIER  MODULUS 

It  is  well  known  that  for  the  one-dimensional  case  the  Fourier  modulus 
is  usually  ambiguous  [11].  That  is,  there  may  be  many  objects  that  have  the 
same  Fourier  modulus.  Going  further,  there  may  even  be  many  real,  nonnega¬ 
tive  objects  that  have  the  same  Fourier  modulus.  One  therefore  might  ques¬ 
tion  why  one  should  even  attempt  to  construct  an  image  from  Fourier  modulus 
data  if  one  could  get  many  different  solutions.  Fortunately  the  one¬ 
dimensional  analysis  does  not  carry  over  into  two  dimensions,  and,  as  will 
be  shown  later,  the  solution  almost  always  is  unique  for  the  two-dimensional 
case.  In  Section  4.1  the  previously  known  theory  is  briefly  reviewed  and 
new  experimental  construction  results  are  shown  in  Section  4.2  for  a  variety 
of  objects. 

4.1  UNIOUENESS  THEORY 

We  first  consider  the  one-dimensional  case,  although  the  problem  in  which 
we  are  most  interested  is  two-dimensional.  Switching  to  continuous  vari¬ 
ables,  the  Fourier  transform  relationship,  of  which  Eq.  (1)  is  a  discrete 
form,  is 

F(u)  -  f  f(x)e_i2lTUX  dx  (12) 

—  00 

If  f(x)  is  of  finite  extent,  that  is,  if  it  is  zero  outside  a  finite  in¬ 
terval  [a,  b],  then  the  function  F(z)  is  an  entire  function  (analytic  every¬ 
where)  of  exponential  type.  The  variable  z  is  a  complex  variable  of  which  u 
is  the  real  part.  Such  functions  are  completely  specified  by  their  complex 
zeroes  z^,  i  «  1,  2,  ...,  satisfying  F(z^)  *  0.  In  particular,  from  the 
Hadamard  factorization  theorem  we  have  [12] 

Ff z)  .  Azq  5  (1  -  z/z.)  (13) 

i-1  1 

where  A  is  a  constant  and  q  is  the  order  of  the  zero  at  z  ■  0. 


Th  ea son  that  the  one-dimensional  case  is  usually  not  unique  is  as 
follow-..  In  theory,  since  F(z)  is  analytic,  given  any  finite  segment  of 
F(z),  F(z)  can  be  determined  everywhere  by  analytic  continuation.  Similarly, 
since  Ffz)  is  analytic,  so  is  F(zlF(z),  which  is  equal  to  I  F( u ) I  ^  on  the 
real  line.  (The  overbar  indicates  complex  conjugation.)  Therefore,  given 
|F(u)  r  it  is  theoretically  possible  to  determine  the  complex  zeroes  of 
F(z)F(F).  If  the  complex  zeroes  of  F(z)  were  known,  then  Eq.  (13)  could  be 
used  to  determine  F(z)  and  in  particular  F(u).  However,  for  each  zero  zi 
of  F(z),  F(z)F(z)  *  0  both  at  z.  and  at  1...  Therefore  for  each  non- 
real-valued  complex  zero  of  F(z)  there  is  a  two-fold  ambiguity  of  the  solu¬ 
tion.  If  there  are  m  complex  zeroes  then  there  would  be  2m-^  different 
solutions,  and  in  general  m  is  infinite.  In  general  it  is  not  known,  how¬ 
ever,  what  number  of  solutions  are  nonnegative. 

Although  the  uniqueness  problem  is  very  severe  in  the  one-dimensional 
case,  it  fortunately  does  not  appear  to  be  a  problem  in  the  two-dimensional 
case.  Bruck  and  Sodir  [13]  analyzed  the  case  of  an  object  made  up  of  an  ar¬ 
ray  of  delta-functions.  The  Fourier  transform  then  reduces  to  a  discrete 
sum,  and  the  discrete  sum  can  be  shown  to  be  equivalent  to  a  polynomial  of 
complex  variables.  The  zeroes  (roots)  of  the  polynomial  are  associated  with 
the  complex  zeroes  discussed  above.  It  can  then  be  shown  that  the  ambiguity 
of  the  solution  is  determined  by  the  number  of  factors  into  which  the  poly¬ 
nomial  can  be  factored.  In  the  one-dimensional  case  a  polynomial  can  always 
be  expressed  as  the  product  of  m  prime  factors,  where  m  is  equal  to  the  order 

of  the  polynomial,  which  is  equal  to  the  number  of  discrete  points  across 

] 

the  object.  Therefore,  for  the  one-dimensional  case  there  are  2  solu¬ 
tions.  On  the  other  hand,  two-dimensional  polynomials  are  rarely  factorable. 
Therefore,  for  the  two-dimensional  case  the  solution  is  usually  unique  (to 
within  a  180*  rotation  and  a  translation). 

It  can  be  argued  that  this  analysis  can  be  extended  to  the  case  of  a  con¬ 
tinuous,  extended  object  (as  opposed  to  a  sampled  object)  by  noting  that  F(z) 
for  a  continuous  object  can  be  closely  approximated  by  a  polynomial  [14], 
Therefore,  although  further  analytical  work  is  required  to  definitively  an¬ 
swer  the  question  of  uniqueness  for  continuous  objects  in  two-dimensions, 
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one  would  expect  the  solution  to  be  generally  unique  (except  for  ver  ,  special 

« 

cases)  in  two  dimensions. 

Fried  [15]  has  considered  the  special  case  for  which  the  Fourier  trans¬ 
form  is  separable  in  its  two  orthogonal  coordinates,  that  is,  it  can  be 
written  as  the  product  of  two  one-dimensional  functions  (and  if  the  Fourier 
transform  is  separable,  then  so  is  the  object  function).  This  case  can  be 
shown  to  have  a  high  degree  of  non-uniqueness.  The  high  degree  of  ambiguity 
is  expected  since  the  separability  condition  changes  the  two-dimensional 
problem  into  two  one-dimensional  problems,  and  the  high  degree  of  ambiguity 
is  well  known  in  one  dimension,  as  discussed  earlier. 


Fried  [15]  goes  one  step  further  and  shows  how  to  generate  an  arbitrary 
number  of  nonnegative  (one-dimensional  or  separable  two-dimensional)  func¬ 
tions  having  the  same  Fourier  modulus.  He  has  since  found  that  the  functions 
generated  by  his  method  tend  to  be  smooth 


4.2  UNIQUENESS  EXPERIMENTS 

The  uniqueness  of  the  solution  was  empirically  tested  on  several  differ¬ 
ent  objects  in  order  to  provide  a  practical  answer  to  the  uniqueness 
question . 

Digitized  photographs  of  several  different  satellites  were  provided  to 
ERIM  by  B.L.  McGlamery  (Visibility  Laboratory,  U.C.  San  Diego).  Since  it 
was  suspected  that  the  uniqueness  property  may  depend  on  the  complexity  of 
the  object,  we  computed  several  versions  of  each  satellite,  each  version  at 
a  different  resolution  (at  a  different  magnification).  The  lower  resolution 
versions  provided  examples  of  less  complicated  objects. 

The  method  of  generating  lower  resolution  versions  was  to  resample  the 
fine  resolution  images  on  a  coarser  grid.  Points  on  the  coarser  grid  that 
fell  between  the  points  on  the  finer  grid  were  given  values  according  to  a 
two-dimensional  separable  triangular  filter  operating  on  the  four  nearest 
neighbor  points  on  the  finer  grid.  This  type  of  resampling  does  not  give  an 
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derragm location  •  it  has  a  tendency  to  overemphasize  some  features 
t  >  object  :>nd  Tins  otr.ers.  However,  for  the  purpose  of  these  experiments 
ar  accurate  demaqnif  i  rati  or.  was  not  reauired;  it  was  only  desired  that  a 

smaller,  less  complicated  object  be  synthesized. 

The  objects  used  are  listed  in  Table  1.  In  order  to  avoid  the  possi¬ 
bility  of  aliasing,  the  objects  were  initially  all  scaled  to  fit  within  a 

square  of  size  64  by  64  pixels,  which  was  imbedded  in  an  array  of  size  128 

by  128.  The  magnification  factors  required  for  the  initial  scaling  are 

listed  in  Table  1.  If  a  telescope  of  diameter  1.6  meters  and  a  wavelength 

of  550  nm  were  used,  and  the  image  sample  size  were  equivalent  to  (x/20)  = 

_7 

550  nm/3.2  m  =  1.719  x  10  radians  *  0.03545  are-sec  per  pixel,  then  the 
apparent  distances  corresponding  to  the  objects  of  size  64  pixels  are  given 
by  (image  scale  in  meters  per  pixel)  ( 2D / x ) ,  which  are  listed  in  tne  last 

column  of  Table  1.  In  addition  to  objects  of  size  64  pixels,  objects  of  size 
32,  16,  8  and  4  pixels  were  generated  using  the  resampling  method  described 

above.  That  is,  a  total  of  40  test  objects  (8  objects  at  5  magnifications 

each)  were  used  in  the  experiments.  The  smaller,  lower-resolution  objects 
can  be  thought  of  as  being  at  proportionally  greater  distances  from  the  tele¬ 
scope.  Because  of  the  way  that  the  resampling  was  done,  the  objects  "of  size 
4  pixels"  are  in  fact  only  3  pixels  across;  in  addition,  except  for  the  Film 
No.  508  and  608  cases,  one  pixel  dominates  the  others;  consequently,  the  ob¬ 
jects  of  size  4  pixels  tend  to  look  very  much  like  single  point  scatterers. 

The  Fourier  modulus  data  for  each  test  object  was  generated  by  Fourier 
transforming  the  test  object,  then  multiplying  the  complex  Fourier  transform 
by  W,  the  weighting  function  due  to  a  circular  telescope  aperture,  and  then 
taking  the  modulus  (or  magnitude).  The  weighting  function  is  given  by  the 
autocorrelation  of  a  circular  aperture  of  diameter  62  pixels  in  the  Fourier 
domain.  For  this  set  of  experiments  the  effects  of  diffraction  due  to  the 
telescope  aperture  were  included,  but  the  Fourier  modulus  data  was  otherwise 
free  of  error.  For  comparison  with  the  construction  results,  the  inverse 
Fourier  transform  of  the  weighted  complex  Fourier  transform  was  computed. 
The  resulting  image,  which  is  referred  to  as  the  reference  object,  is  equiv¬ 
alent  to  a  diffraction-limited  image  of  the  object. 
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The  same  set  of  340  iterations  as  discussed  in  Section  2  was  us  each 

of  the  forty  test  cases.  Figures  10  to  17  show  the  reference  obje.  ,  and  a 
corresponding  constructed  image  for  all  forty  cases.  Images  that  came  out 
inverted  were  re-inverted  in  order  to  take  on  the  same  orientation  as  the 
reference  object. 

A  problem  that  occurred  was  that  the  algorithm  usually  did  not  converge 
completely  to  a  solution:  it  usually  stagnated  at  a  local  minimum  at  which 
it  produced  an  image  of  qood  quality  but  with  a  set  of  stripes  superimposed 
across  the  image.  When  the  algorithm  stagnated  at  such  a  local  minimum,  it 
was  not  possible  to  determine  with  certainty  whether  the  solution  it  was  con¬ 
verging  toward  was  exactly  the  same  as  the  reference  object.  As  discussec 
in  Section  3,  the  stagnation  at  an  output  image  having  stripes  is  most  likely 
a  problem  of  algorithm  stagnation  at  a  local  minimum  of  EQ  and  does  not 
represent  a  uniqueness  problem.  In  almost  all  of  the  different  constructed 
images  shown  in  Figures  10-17,  the  primary  difference  between  the  constructed 
image  and  the  reference  object  is  the  presence  of  the  stripes.  Therefore  we 
conclude  from  those  cases  that  if  the  algorithm  were  improved  in  such  a  way 
as  to  avoid  stagnation  with  stripes,  then  it  would  probably  converge  to  ar. 
output  image  identical  to  the  reference  object.  That  is,  these  results  imply 
that  for  most  cases  the  solution  is  unique. 

The  stripes  are  usually  more  difficult  to  see  when  reproduced  than  they 
are  when  the  image  is  viewed  on  a  T.V.  monitor,  which  has  greater  dynamic 
range  than  a  paper  print  of  the  image.  Usually  there  is  a  single  set  of 
stripes  across  the  image.  In  any  given  image  the  stripes  have  an  average 
period  and  angular  orientation,  but  across  the  image  the  stripes  vary  signif¬ 
icantly  about  those  averages.  In  some  cases,  more  than  one  set  of  stripes 
appear  in  a  given  image;  for  example,  in  Fiaure  12(B)  there  are  strong  (high 
contrast)  stripes  both  horizontally  and  at  about  a  45  degree  angle.  The 
stripes  extend  over  the  entire  f ield-of-view  but  have  maximum  contrast  over 
the  extent  of  the  object.  The  fact  that  the  stripes  extend  beyond  the  known 
extent  of  the  object  (which  is  half  the  diameter  of  the  autocorrelation)  will 
probably  enable  us,  some  time  in  the  future,  to  develop  an  algorithm  to  de¬ 
tect  them  and  eliminate  them. 


38 


'  lists  the  rms  error  EQ,  as  defined  by  Eq.  (5),  for  the  forty 
Co.-;?  t-  .. ted  images.  This  is  a  measure  of  how  well  the  constructed  image 
agrees  with  the  simultaneous  constraints  of  nonneaati vity  and  having  a 
Fourier  modulus  equal  to  the  Fourier  modulus  of  the  reference  object;  that 
is,  it  is  a  measure  of  how  close  it  is  to  a  solution.  For  a  given  size  of 
object,  constructed  images  with  larger  values  of  EQ  are  known  to  be 
farther  from  a  solution.  Therefore,  for  example,  the  32-pixel  size  and 
16-pixel  size  constructed  images  of  Film  No.  302  shown  in  Figures  12(0)  and 
(F)  should  probably  be  ignored  because  they  are  much  farther  from  solutions 
than  are  the  constructed  images  for  the  other  objects  (and  so  they  are  not 
necessarily  an  indication  of  non-uniqueness). 

If  the  cases  that  did  not  converge  very  close  to  a  solution  (having  large 
values  of  EQ)  are  ignored,  and  if  perturbations  due  to  stripes  are  ig¬ 
nored,  then  most,  but  not  all,  of  the  image  construction  results  shown  in 
Figures  10-17  seem  to  indicate  that  the  solutions  are  unique.  The  most  dis¬ 
turbing  departure  from  uniqueness  came  in  the  case  of  the  64-pixel  size  of 
Film  No.  401,  FLTSATCOM.  The  constructed  image  has  the  same  overall  shape 
as  the  reference  object  but  differs  from  it  in  its  details,  despite  the  fact 
that  for  this  object  an  additional  380  iterations  were  performed.  Subsequent 
runs  of  the  algorithm  using  different  starting  inputs  encountered  similar 
problems  for  this  object.  From  Figure  14(A)  one  can  see  that  for  this  par¬ 
ticular  object  its  bright  central  part  is  nearly  circularly  symmetric.  Since 
a  circularly  symmetric  object  can  be  described  as  a  function  of  a  single 
variable,  its  radial  coordinate,  one  could  hypothesize  that  such  objects 
would  have  the  same  high  degree  of  ambiguity  as  one-dimensional  objects,  but 
this  is  not  known.  Therefore,  although  most  objects  of  interest  are  probably 
unique,  there  appear  to  be  special  classes  of  objects  for  which  the  solution 
is  not  unique.  From  the  results  shown  in  Figures  10-17  it  appears  that  the 
non-unique  (ambiauous)  cases  are  the  exceptions  rather  than  the  rule. 


TABLE  1 

OBJECTS  FOR  UNIQUENESS  STUDY 


F  ilm 
No. 

Object  Description 

Unsealed 

Pixel 

Size 

Magn if ication 
for  64-Pixel 
Width 

Apparent  Distance 
for  64-Pixel 

Width  (km) 

204 

P72-2  Sensor  Testbed  (A) 

0.0174m 

0.267 

379 

205 

P72-?  Sensor  Testbed  (B) 

0.0174m 

0.267 

379 

302 

777  Communications  Sat(A) 

0.0174m 

0.267 

379 

314 

777  Communications  Sat(B) 

0.0174m 

0.279 

363 

401 

FLTSATCOM 

0.0581m 

0.250 

1352 

509 

LST  (A) 

0.0871m 

0.356 

1423 

508 

LST  (B) 

0.0871m 

0.291 

1741 

608 

SURVSATCOM 

0.116m 

0.356 

1896 

TABLE  2 

RMS  ERROR  E0  OF  CONSTRUCTION  RESULTS  FOR 
UNIQUENESS  STUDY 


Object 

Object 

.  Size  (Pixels) 

Film  No.) 

64 

32 

16 

8 

4 

204 

0.0327 

0.0234 

0.00914 

0.00051 

0.00048 

205 

0.0262 

0.00421 

0.00074 

0.00058 

0.00049 

302 

0.0189 

0.0553 

0.0154 

0.00071 

0.00053 

314 

0.0156 

0.00268 

0.00813 

0.00164 

0.00055 

401 

0.00240 

0.00101 

0.00061 

0.00050 

0.00048 

509 

0.00782 

0.00660 

0.00072 

0.00053 

0.00048 

508 

0.00296 

0.00179 

0.00197 

0.00071 

0.00048 

608 

0.00368 

0.00169 

0.00073 

0.00100 

0.00054 
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Figure  10.  Reference  objects  and  ccnslr in, ages  for  ^ .  I No.  204, 
P72-2(A>.  reference  object:  and  construe.;.':,  -rage  for-  f  /:-4-pixe1 
width,  (C,  0)  3?- pixel  w  i  dt  h ,  'c,  m  Ifc-pixel  width,  !r:y  ?1(  d-pixel 

wid4h,  ana  (I,  J;  4-p'i^e!  widtv. 
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Figure  11.  Reference  objects  and  constructed  images  for  Film  No.  205, 
P72-2(B).  Reference  object  and  constructed  image  for  (A,  B)  64-pixel 
width,  (C,  0)  32-pixel  width,  (E,  F)  16-pixel  width,  (G,  H)  8-pixel 

width,  and  (I,  J)  4-pixel  width. 
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Figure  12.  Reference  objects  and  constructed  images  for  Film  No.  302, 
777(A).  Reference  object  and  constructed  Image  for  (A,  B)  64-pixel 
width,  (C,  D)  32-pixel  width,  (E,  F)  16-pixel  width,  (G,  H)  8-pixel 

width,  and  (I,  J)  4-pixel  width. 
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Figure  13.  Reference  objects  and  constructed  images  for  Film  No.  314, 
777(B).  Reference  object  and  constructed  image  for  (A,  B)  64-pixel 
width,  (C,  D)  32-pixel  width,  (E,  F)  16-pixel  width,  (G,  H)  8-pixel 

width,  and  (I,  J)  4-pixel  width. 
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figure  14.  Reference  objects,  and  a>n:.Vw.ted  i  for  F  i ;  :v, 

FLTSATCOM.  Reference  object  rind  cr-v. trusted  i  c  ■  :  {,\.  b 

width,  r,  0)  32-pixel  width,  if.  ;6-pi,.‘:l  wo.--,  i) 

w  /  J  t  h»  d  r\  d  '  ] ,  . 1 '  t  - 1)  7  x  el  w  i  J  r  v  . 


Figure  15.  Reference  objects  and  constructed  images  for  Film  No.  509, 
LST(A).  Reference  object  and  constructed  image  for  (A,  B)  64-pixel 
width,  (C,  D)  32-pixel  width,  (E,  F)  16-pixel  width,  (G,  H)  B-pixel 

width,and  (I,  J)  4-pixel  width. 
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Figure  16.  Reference  objects  and  constructed  images  for  Film  No.  508, 
LST(B).  Reference  object  and  constructed  image  tor  (A,  B)  64-pixel 
width,  (C,  D)  32-pixel  width,  (E,  F)  16-pixel  width,  (G,  H)  8-pixel 

width, and  (I,  J)  4-pixel  width. 


Figure  17.  Reference  objects  and  constructed  images  for  Film  No.  608, 
SURVSATCOM.  Reference  object  and  constructed  image  for  (A,  B)  64-pixel 
width,  (C,  D)  32-pixel  width,  (E,  F)  16-pixel  width,  (G,  H)  8-pixel 

width, and  (I,  J)  4-pixel  width. 
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ALGORITHM  CONVERGENCE 


Several  different  aspects  of  the  convergence  of  the  algorithm  are  of  in¬ 
terest.  A  few  variations  of  the  algorithm  were  compared  for  convergence 
speed.  Investigations  were  made  into  the  convergence  speed  as  a  function  of 
object  complexity  and  of  Fourier  modulus  signal-to-noise  ratio.  Attempts 
were  made  to  cause  the  solution  to  move  away  from  a  local  minimum  character¬ 
ized  by  the  stripes  discussed  in  Sections  3  and  4. 

5.1  COMPARISON  OF  DIFFERFNT  ALGORITHMS 

Four  different  versions  of  the  algorithm  were  compared.  All  four  ver¬ 
sions  use  the  same  first  three  steps  described  in  Section  2.  They  differ  in 
the  method  in  which  the  new  input  gK+^  is  chosen.  The  first  method,  which 
we  call  the  error-reduction  approach,  is  given  by  Eq.  (2)  in  Section  2: 

|9k,mn'  <"•  ">*  •> 

(I4J 

0,  (m,  n)  e  y 

where  g^  mn  and  g'k  mn  are  the  input  and  output,  respectively,  at  the 

kth  iteration  and  y  is  the  set  of  points  at  which  q\  violates  the 

3  k ,  mn 

object-domain  constraints.  The  principal  object-domain  constraint  is  that 
g'  be  nonnegative  everywhere;  an  auxiliary  constraint  is  that  the  extent 
(diameter)  of  the  object  not  exceed  the  known  diameter  which  is  half  the 

diameter  of  the  given  autocorrelation.  For  the  error-reduction  approach  the 

new  input  is  equal  to  the  current  output  modified  to  satisfy  the  object- 
domain  constraints.  Satisfying  the  object  domain  constraints  consists  of 
setting  the  output  to  zero  wherever  it  violates  the  constraints.  In  Appendix 
A  is  a  proof  that  the  error  reduction  approach  converges  in  the  sense  that 
the  mean  squared  error,  as  defined  in  Eqs.  (4)  and  (5)  in  Section  2,  is  mono- 
tonically  decreasing  with  each  successive  iteration. 

In  addition  to  the  error-reduction  approach,  three  different  versions  of 
the  more  general  input-output  approach  were  investigated.  The  basic  idea 
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behind  the  input-output  approach  is  that  the  first  three  operations  of  the 
algorithm  —  Fourier  transformation,  satisfaction  of  Four ier-domain  con¬ 
straints,  and  inverse  Fourier  transformation  —  can  be  viewed  as  a  nonlinear 
system  having  an  input  g  and  an  output  g'.  By  appropriate  changes  in  the 
input  one  attempts  to  drive  the  output  to  satisfy  the  object-domain  con¬ 
straints.  Since  the  output  by  definition  satisfies  the  Four ier-domain  con¬ 
straints,  it  is  a  solution  to  our  problem  if  it  can  be  made  to  satisfy  the 
object-domain  constraints. 


For  a  problem  very  similar  to  the  one  at  hand,  it  has  been  shown  that 
for  small  changes  in  the  input,  the  expected  value  of  the  corresponding 
change  in  the  output  is  a  constant,  a,  times  the  change  in  the  input  [4]. 
Since  additional  nonlinear  terms  also  appear  in  the  output,  the  change  in 
the  output  due  to  a  particular  change  in  the  input  cannot  be  predicted  ex¬ 
actly.  Nevertheless,  by  appropriate  changes  in  the  input,  the  output  can  be 
pushed  in  the  general  direction  that  is  desired.  If  a  change  Agmn  is  de¬ 
sired  in  the  output,  then  a  logical  choice  for  a  change  in  the  input  to 
achieve  that  change  in  the  output  would  be  6Agmn,  where  8  is  a  constant 
ideally  equal  to  a~  .  For  the  problem  at  hand,  the  desired  change  in  the 
output  is 

0,  (m,  n) $  y 

Ag  =  (15) 

3mn 

-q'  ,  (m,  n)  e  y 

'  ymn’  '  ’  '  ' 


that  is,  where  the  constraints  are  satisfied  one  does  not  require  a  change 
of  the  output,  but  where  the  constraints  are  violated,  the  desired  change  in 
the  output  is  one  that  drives  it  to  a  value  of  zero  (and  therefore  the  de¬ 
sired  change  is  the  negative  of  the  output  at  those  points).  Therefore  a 
logical  choice  for  the  next  input  is  g  +  BAg,  that  is, 

I  9k, nr-  (m-n)<t  ’ 

9k+l,im  = 

«k, mn  -  89'k,„,. <■.">  «  Y 


We  will  refer  to  the  use  of  Eq.  (16)  as  the  basic  input-output  approach. 
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An  interesting  property  of  the  nonlinear  system  (consisting  of  the  three 
steps  mentioned  above)  is  that  if  an  output  g'  is  used  as  an  input,  then  its 
output  will  be  itself.  Since  the  Fourier  transform  of  g'  already  satisfies 
the  Four ier-domain  constraints,  g'  is  unaffected  as  it  goes  through  the  sys¬ 
tem.  Therefore,  irrespective  of  what  input  actually  resulted  in  the  output 
g ' ,  the  output  g1  can  be  considered  to  have  resulted  from  itself  as  an  input. 
From  this  point  of  view,  another  loaical  choice  of  the  next  input  is 

Sk.mn*  <*•">*» 

9k+l,nm  *  n?) 

9k, mn  '  fl9'k,mn’  *  * 

We  will  refer  to  the  use  of  Fq.  (17)  as  the  output-output  approach. 

Note  that  if  e  ■  1  in  Eq.  (17),  then  the  output-output  approach  reduces 
to  the  error-reduction  approach  of  Eq.  (14).  Since  the  optimum  value  of  e 
Is  usually  not  unity,  the  error-reduction  approach  can  be  looked  on  as  a  sub- 
optimal  version  of  a  more  general  approach. 


A  fourth  method  of  choosing  the  next  input  which  we  investigated  is  a 
combination  of  the  upper  line  of  Fo.  M7)  with  the  lower  line  of  Fq.  M6): 


I  , mn ’ 


Vi, 


mn 


(m,n)  i 


9k,mn  -  s9'k,mn*  <•.">  «  * 
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We  will  refer  to  the  use  of  Fq.  (18)  as  the  hybrid  input-output  approach. 
The  hybrid  input-output  approach  is  an  attempt  to  avoid  a  stagnation  problem 
that  tends  to  occur  with  the  output-output  approach.  The  output-output  ap¬ 
proach  often  works  itself  into  a  situation  in  which  the  output  on  successive 
iterations  does  not  change,  despite  being  far  from  a  solution.  For  the  hy¬ 
brid  input-output  approach,  on  the  other  hand,  if  at  a  given  pixel  (m,  n) 
the  output  remains  neoative  for  more  than  one  iteration,  then  the  corre¬ 
sponding  point  in  the  input  continues  to  grow  larger  and  larger  until  even¬ 
tually  that  output  point  must  go  nonnegative. 
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The  four  approaches  discussed  above  were  compared  by  using  them  on  the 
same  Fourier  modulus  data  and  with  the  same  starting  input.  For  each  ap¬ 
proach  several  different  values  of  the  parameter  e  were  tried.  The  principal 
problem  with  the  error-reduction  approach  is  that  it  tends  to  stagnate  after 
a  few  iterations.  The  starting  point  for  the  iterations  was  chosen  to  be  a 
partially  constructed  image  of  the  Film  No.  509,  LST(A)  object  (at  a  64-pixel 
size),  on  which  the  error-reduction  approach  was  making  very  slow  progress. 
Ten  iterations  of  each  of  the  four  approaches  followed  by  ten  iterations  of 
the  error-reduction  approach  (that  is,  a  total  of  twenty  iterations)  were 
performed  using  that  same  starting  input.  The  reason  that  each  approach  was 
followed  by  ten  iterations  of  the  errorreduction  approach  is  that  in  some 
cases  definite  progress  is  being  made  with  an  input-output  approach  even 
though  the  rms  error  gets  worse  with  each  iteration.  The  relationship  be¬ 
tween  the  rms  error  and  the  visual  image  quality  is  not  fully  understood, 
although  of  course  one  would  expect  a  high  degree  of  correlation  between  the 
two.  For  the  approaches  for  which  the  rms  error  does  not  seem  to  be  a  good 
indication  of  the  image  quality,  we  found  that  the  rms  error  could  be  made 
to  be  a  good  measure  of  the  image  quality  by  performing  a  few  (say  ten)  iter¬ 
ations  of  the  error-reduction  approach  on  the  results  of  the  input-output 
approach . 

Figure  18  shows  a  plot  of  the  rms  error  after  the  twenty  iterations  for 
each  of  the  input-output  approaches  as  a  function  of  the  parameter  b.  Recall 
that  the  output-output  approach  with  e  =  1.0  is  equivalent  to  the  error 
reduction  approach.  Figure  18  shows  that  the  hybrid  input-output  approach 
is  superior  to  the  others  in  this  case,  and  that  the  optimal  value  of  B  is 
about  unity. 

The  manner  in  which  the  diameter  constraint  (limiting  the  diameter  of 
the  constructed  image  to  half  the  diameter  of  the  autocorrelation)  was  im¬ 
posed  was  found  to  have  a  sianificant  impact  on  the  convergence  speed.  At 
the  beginning  of  this  program  the  strategy  for  applying  the  diameter  con¬ 
straint  was  to  impose  it  very  loosely  ( i . e . ,  allow  the  image  to  have  a  larger 
diameter)  for  the  early  iterations,  then  tighten  up  the  diameter  for  later 
iterations  after  the  image  distribution  became  better  confined  to  an  area  of 
the  desired  size.  This  strategy  helped  to  avoid  the  previously  encountered 
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Algorithm  Parameter  fi 


Figure  18.  Rms  error,  £Q,  after  sequence  of  20  Iterations  versus 

algorithm  parameter  e  for  the  basic  Input-output  approach  (a),  the 
output-output  approach  (•)  and  the  hybrid  Input-output  approach  (x) 


problem  of  having  the  diameter  constraint  inadvertently  chopping  off  an  edge 
of  the  desired  image.  It  was  discovered  that  reversing  the  sequence  resulted 
in  an  increase  in  convergence  speed:  apply  a  very  tight  diameter  constraint 
for  the  early  iterations,  and  loosen  that  constraint  as  the  iterations  pro¬ 
gress.  This  strategy  helps  to  move  energy  from  the  regions  clearly  outside 
the  object's  support  to  within  the  objects  support  during  the  early  itera¬ 
tions.  This  results  in  a  convergence  speed-up  by  as  much  as  a  factor  of  two. 

5.2  DATA  DEPENDENCY  OF  CONVERGENCE 

The  convergence  speed  of  the  algorithm  depends  both  on  the  complexity  of 
the  object  and  upon  the  signal-to-noise  ratio  of  the  Fourier  modulus  data. 
The  most  convenient  scale  by  which  convergence  is  measured  is  the  rms  error 
Eq  or  Ep,  both  of  which  are  consistently  meaningful  only  after  a  few 

iterations  of  the  error  reduction  approach  are  performed. 

The  data  generated  during  the  studies  of  the  noise  properties  of  the  con¬ 
structed  images  show  that  the  number  of  iterations  required  was  the  least 
for  the  data  of  the  lowest  signal-to-noise  ratio.  For  the  data  having  the 
very  worst  signal-to-noise  ratio  (photoelectron  scale  factor  =  6),  only  a 
few  iterations  were  required  for  convergence.  That  is  to  say,  after  just  a 
few  iterations  no  further  progress  could  be  made  although  EQ  was  still 
large,  and  the  solution  had  gone  as  far  as  it  usefully  could.  For  the  higher 
signal-to-noise  ratio  cases,  additional  iterations  further  reduced  the  rms 
error,  so  in  that  sense  the  convergence  time  was  longer. 

The  data  generated  durinq  the  studies  of  the  uniqueness  of  the  solution 
showed  that  the  objects  of  lesser  complexity  (smaller  magnification)  con¬ 
verged  to  lower  values  of  rms  error  E  after  fewer  iterations  than  the 

0 

more  complex  objects.  However,  the  total  number  of  iterations  performed  be¬ 
fore  the  algorithm  stagnated  showed  only  a  weak  dependency  on  object  complex¬ 
ity.  For  the  cases  of  low  object  complexity  stagnation  occurred  at  consid¬ 
erably  lower  values  of  EQ  than  for  cases  of  higher  object  complexity,  but 

stagnation  tended  to  occur  sooner  for  cases  of  low  object  complexity. 
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b.3  THE  STRIPES  phenomenon 


As  mentioned  in  Sections  3  and  4  the  algorithm  often  stagnates  at  a  local 
minimum  characterized  by  a  pattern  of  stripes  across  the  image.  In  most 
cases  the  stripes  are  barely  noticeable  and  are  of  low  contrast  superimposed 
on  an  otherwise  excellent  constructed  image.  In  some  cases  the  stripes  are 
of  high  enough  contrast  to  be  objectionable,  although  they  still  permit  the 
object  to  be  recognized.  As  shown  in  Section  3,  different  starting  inputs 
converqe  to  images  having  stripes  of  different  angular  orientations  and 
spatial  frequencies  (periods).  There  seems  to  be  a  tendency  for  the  stripes 
to  be  oriented  at  the  same  angle  as  a  prominent  feature  of  the  object.  The 
stripes  extend  over  the  entire  f ield-of-view,  beyond  the  known  diameter  of 
the  object.  The  contrast  of  the  stripes  tends  to  be  greater  over  the  extent 
of  the  object  and  less  outside  the  extent  of  the  object,  The  stripes  appear 
only  when  the  signal-to-noise  ratio  is  high. 

We  believe  that  the  existence  of  the  stripes  does  not  represent  a  funda¬ 
mental  uniqueness  problem.  When  the  algorithm  stagnates  on  an  output  having 
stripes,  the  rms  errors  FQ  and  Ep  are  never  zero,  and  so  it  is  known 
that  it  is  not  exactly  a  solution,  although  it  is  usually  close  to  a  solu¬ 
tion.  Tf  the  algorithm  happens  to  find  the  solution  free  of  stripes,  then 
the  rms  error  EQ  is  lower  than  for  the  case  having  stripes,  providing  a 
way  to  distinguish  the  correct  solution  from  the  rolution  having  stripes. 
In  addition,  since  a  valid  solution  cannot  have  a  diameter  exceeding  half 
the  diameter  of  the  autocorrelation,  the  fact  that  the  stripes  extend  beyond 
that  diameter  provides  a  way  to  distinguish  the  correct  solution  from  a  so¬ 
lution  having  stripes. 

For  the  one-dimensional  case  it  is  known  that  if  a  single  pair  of  complex 
zeroes  is  moved  (to  locations  other  than  their  conjugates)  around  symmetri¬ 
cally  in  the  Fourier  complex  plane,  then  the  result  in  the  image  plane  can 
be  a  function  similar  to  the  original  but  having  a  sinusoid  added  over  the 
interval  of  its  support  fl2l.  This  particular  operation  results  in  a  Fourier 
transform  with  a  Fourier  modulus  different  from  the  original  Fourier  modulus; 
therefore  the  new  object  created  does  not  cause  a  problem  of  uniqueness. 
Furthermore,  this  one-dimensional  analysis  does  not  carry  over  directly  into 
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two  dimensions.  Nevertheless,  this  one-dimensional  phenomenon  might  ulti¬ 
mately  be  connected  in  some  way  with  the  stripes  phenomenon. 

An  attempt  was  made  at  moving  away  from  the  constructed  image  having 
stripes.  After  the  algorithm  converged  to  an  image  having  stripes  and  seemed 
to  stagnate  there,  random  noise  was  added  to  the  image  and  then  further  iter¬ 
ations  were  performed.  It  was  found  that  the  algorithm  quickly  converged 
back  to  an  image  with  essentially  the  same  stripes,  even  when  large  amounts 
of  noise  had  been  added.  Only  by  starting  with  a  different  initial  input 
having  no  stripes  present  could  one  converge  to  an  image  having  a  different 
set  of  stripes.  Further  work  will  be  necessary  to  develop  a  method  of  break¬ 
ing  away  from  stagnation  at  local  minima  having  stripes. 
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6 

CONCLUSIONS  AND  RECOMMENDATIONS  FOR  FUTURE  RESEARCH 


The  results  shown  here  indicate  that  construction  of  images  by  the  iter¬ 
ative  method  from  Fourier  modulus  data  provided  by  stellar  speckle  interfer¬ 
ometry  should  prove  to  be  successful  for  imaging  earth-orbiting  satellites. 
With  the  levels  of  photon  noise  present  in  practical  situations  it  should  be 
possible  to  construct  images  that  are  diffraction-limited  in  resolution. 
The  solution  can  generally  be  counted  on  to  be  unambiguous  (i.e.,  unique) 
except  for  special  classes  of  objects,  for  example,  objects  that  can  be 
described  as  separable  functions.  The  image  construction  algorithm  has  been 

developed  to  a  point  were  convergence  requires  about  100  iterations,  which 

takes  less  than  two  minutes  on  an  array  processor  (Floating  Point  System 
AP-1P0B)  for  array  sizes  128  by  128. 

More  particular ly,  the  major  results  of  this  research  effort  are  as  fol¬ 
lows.  The  rms  error  of  the  constructed  images  increases  in  a  gradual  and 
predictable  manner  as  the  rms  error  of  the  Fourier  modulus  estimate  in¬ 

creases,  roughly  as  the  square  root  of  the  rms  error  of  the  Fourier  modulus 
estimate.  For  photon  noise  levels  up  to  ten  times  worse  than  the  baseline 
case,  the  quality  of  the  constructed  image  is  good.  The  most  recent  theory 
of  the  uniqueness  of  the  solution  suggests,  and  our  experimental  results  ap¬ 
pear  to  confirm,  that  the  solution  is  usually  unique  for  two-dimensional  ob¬ 
jects,  both  for  objects  of  high  complexity  (large  space-bandwidth  product) 
and  low  complexity.  These  results  were  obscured,  however,  by  an  inability 
to  converge  completely  to  a  solution  in  many  cases.  The  algorithm  has  a 

tendency  to  stagnate  at  a  local  minimum  characterized  by  a  set  of  stripes 
across  an  otherwise  recognizable  image  of  the  object.  As  the  signal-to-noise 
ratio  of  the  Fourier  modulus  data  decreases,  the  ambiguity  of  the  solution 
increases,  but  that  ambiguity  takes  the  form  of  noise  in  the  constructed  im¬ 
age  rather  than  a  complete  change  in  the  basic  characteristics  of  the  image. 
Fewer  iterations  of  the  algorithm  were  needed  in  the  low  signal-to-noise 
ratio  cases,  but  the  images  were  correspondingly  poorer.  Images  of  greater 
complexity  tended  to  take  a  larger  number  of  iterations  than  images  of  lower 
complexity.  The  convergence  speech  of  the  algorithm  was  increased  by  the 
manner  in  which  the  diameter  constraint  was  imposed.  A  comparative  study  of 
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a  few  different  versions  of  the  algorithm  showed  that,  among  the  approaches 
tried,  the  hybrid  input-output  approach  converges  the  fastest  by  a  wide 
marqin. 

The  most  critical  questions  have  been  answered  successfully  by  this  re¬ 
search,  but  important  aspects  of  the  problem  remain  to  be  studied. 

Althouoh  useful  imaqes  are  obtained  in  spite  of  it,  the  stripes  phenom¬ 
enon  remains  as  a  puzzling  hindrance  to  the  complete  convergence  of  the  al¬ 
gorithm  to  a  solution.  Theoretical  analysis  is  needed  to  understand  why  the 
stripes  occur  and  how  they  might  be  deali  with.  Further  effort  should  be 
expended  to  determine  ways  of  automatically  identifying  the  presence  of  the 
stripes  and  having  the  algorithm  reject  them.  Once  this  problem  is  solved, 
the  uniqueness  experiments  should  be  repeated  in  order  to  determine  with 
greater  certainty  whether  the  solution  is  usually  unique. 

A  potential  source  of  error  not  addressed  in  this  research  project  is 
the  estimation  of  the  MTF  of  the  speckle  interferometry  process.  Computer 
simulations  should  be  performed  in  order  to  determine  the  effect  of  errors 
in  the  MTF  estimate  on  the  constructed  image  and  whether  they  affect  the 
uniqueness  of  the  solution.  In  addition,  there  may  be  situations  for  which 
it  is  not  always  convenient  or  practical  to  obtain  an  estimate  of  the  speckle 
MTF  from  a  reference  star.  Therefore,  alternative  methods  of  estimating  the 
speckle  MTF  should  be  studied  and  compared,  including  model  fitting,  the 
Worden  subtract  method,  and  the  clipping  method  [10]. 

Further  development  of  the  algorithm  itself  is  needed  to  make  it  converge 
in  real  time.  More  variations  of  the  algorithm  should  be  compared  over  a 
wider  range  of  circumstances.  In  addition,  major  modifications  of  the  algo¬ 
rithm  may  be  called  for.  A  particularly  promising  line  of  attack  would  be 
to  apply  the  discipline  of  control  theory  to  the  input-output  system  concept. 

Image  construction  experiments  should  be  carried  out  on  actual  telescope 
data  gathered  on  orbiting  satellites  to  demonstrate  its  usefulness  in  solving 
the  space  object  identification  problem  and  to  discover  what  additional  prob¬ 
lems  need  to  be  overcome  when  handling  data  from  a  real  sensor. 

Comparison  of  the  imaqe  construction  results  from  noisy  Fourier  modulus 
data  should  be  made  both  with  (computer  simulated)  images  constructed  by  the 
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Knox-Thompson  method  [161  and  with  images  from  a  compensated  imaqinq  system. 
As  the  reflectivity  of  the  object  decreases,  there  is  some  point  beyond  which 
it  is  no  longer  possible  to  accurately  sense  the  wavefront  deformation,  and 
the  compensated  imaging  system  would  no  longer  adequately  correct  for  atmo¬ 
spheric  turbulence.  However,  under  the  same  severely  photon-limited  case, 
speckle  interferometry  can  continue  to  function  since  it  integrates  ovei  many 
frames  and  can  ultimately  have  a  hiqher  siqnal-to-noise  ratio  than  the  com¬ 
pensated  imaging  system  which  must  correct  for  turbulence  based  on  measure¬ 
ments  made  over  an  interval  of  time  of  only  a  few  milliseconds.  Computer 
simulations  should  be  carried  out  on  the  Knox-Thompson  method  using  the  same 
simulated  data  that  was  generated  for  the  noise  properties  study  described 
in  this  report.  Then  using  existing  computer-simulated  results  for  the 
Fourier  modulus  image  construction  method  and  for  compensated  imaqing,  a  com¬ 
parison  of  the  three  methods  should  be  made  for  various  levels  of  photon 
noise.  It  may  also  be  desirable  to  vary  the  turbulence  parameter  r  .  To 
make  the  comparison  complete,  the  same  type  of  Wiener  filtering  should  be 
performed  on  the  Knox-Thompson  and  Fourier  modulus  image  construction  results 
as  is  done  on  the  images  from  the  compensated  imaging  system. 

Methods  of  combining  the  image  construction  from  speckle  interferometry 
Fourier  modulus  data  with  compensated  imaging  systems  should  be  explored. 
It  may  be  possible  to  have  a  compensated  imaqing  system  utilize  the  iterative 
method  for  those  instances  where  low  light  levels  hinder  the  ability  to  track 
the  wavefront  deformation  caused  by  the  turbulent  atmosphere.  Experiments 
both  with  simulated  compensated  images  and  eventually  with  real  data  from  a 
compensated  imaging  system  would  be  needed.  Most  simply,  the  iterative 
method  can  be  thought  of  as  another  post-processing  step  to  obtain  improved 
imagery  from  the  compensated  imaging  system.  When  compensated  images  from 
the  system  being  developed  at  AMOS  on  Maui  become  available,  if  they  are  not 
diffraction-limited  they  should  be  additionally  post-processed  using  the 
iterative  algorithm. 

For  the  case  of  tumbling  or  rotating  objects,  methods  will  have  to  be 
developed  to  take  these  effects  into  consideration.  For  general  time-varying 
objects  It  may  be  possible  to  "lock  on"  to  the  image:  after  acquiring  an  im¬ 
age  at  a  given  time,  as  the  object  changes  and  new  data  is  collected,  only  a 
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few  iterations  of  the  algorithm  may  be  required  to  converge  on  the  new  image 
if  information  about  the  previous  image  is  utilized. 


Methods  of  increasino  the  signal-to-noise  ratio  of  the  data  to  ensure 
good  quality  constructed  images  should  be  developed.  One  method  of  increas¬ 
ing  the  siqnal-to-noise  ratio  is  to  extend  the  spectral  bandwidth  of  speckle 
interferometry.  A  present  limitation  on  the  spectral  bandwidth  is  the  wave¬ 
length  dependence  of  the  impulse  response  for  a  given  optical  path  difference 
caused  by  the  atmosphere.  It  may  be  possible  to  create  a  wavelength- 
independent  impulse  response  using  an  optical  system  composed  of  wavelength- 
sensitive  optics,  such  as  holographic  optical  elements. 

Finally,  the  iterative  method  could  be  applied  to  a  number  of  other  prob¬ 
lems  of  interest  to  the  Air  Force.  It  can  be  applied  to  any  reconstruction 
problem  for  which  only  partial  information  is  available  about  an  object, 
wavefront,  or  signal  and  only  partial  information  is  available  about  its 
Fourier  transform.  Such  applications  include,  in  addition  to  Fourier  modulus 
image  construction,  wavefront  sensing,  spectral  extrapolation,  and  X-ray 
crystal lography  phase  retrieval.  The  iterative  method  can  also  be  applied 
to  synthesis  problems  for  which  one  wants  both  a  function  and  its  Fourier 
transform  to  satisfy  a  given  set  of  constraints  or  have  certain  desirable 
properties.  Such  applications  include  spectrum  shaping  and  the  design  of 
lens  pupil  functions,  antenna  array  phases,  radar  signals,  and  digital 
filters. 
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APPENDIX  A 
CONVERGENCE  PROOF 


In  this  appendix  it  is  shown  that  for  the  error-reduction  approach  the 
error  monotonically  decreases  at  each  iteration, 
t  h 

The  k  iteration  begins  with  gk  m ,  an  estimate  of  the  object  that 
satisfies  the  object-domain  constraints.  First  gk  mn  is  Fourier  trans¬ 
formed,  resulting  in  Gk  pq  =  |  Gk  pq  I  exp^|<  pq)-  Then  in  the  F°urier 
domain,  Gk>pq  is  replaced  by  G*k>pq  =  |Fpq’|  exp(i«k>pq).  The  unnor¬ 
malized  mean-squared  error  is  defined  in  the  Fourier  domain  as 


EFk  *  77  X)(|Sk,pql“  lFpq  P 


(Al) 


pq 

where  p,q  »  0,  1,  ...,  N-l,  which  can  be  expressed  as 


EFk  -  3  2  |Nk.PO 

pq 


(A2) 


where 


Nk,pq  “  Gk,pq  “  Gk,pq  "  ^  lFpq  I  “  I  Gk,pq  P  exp(l^k,pq^ 


(A3) 


Forming  the  new  Fourier  transform  G'k  pq  and  inverse  Fourier  transforming 
results  in 


dk,mn  " ^Gk,pq^  “  pk,mn  '  "k,mn 


+  n, 


(A4) 


where  "k.«n  •^'■1C,lk,pg)- 

The  new  estimate  of  the  object  is  formed  by  setting  g'  equal  to 

k  ^  mn 

zero  wherever  it  violates  the  constraints: 


9k, m-  '"•n>  i  » 


sk»1. 


mn 


(A5) 


0,  (m,n)  e  y 
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Let  g.  ..  be  expressed  in  terms  of  g,,  „  as  follows 
^k+l ,mn  K,mn 


^k+l ,mn  “  ^k,mn  +  dk,mn  “  ®k,mn  nk,mn  dk,mn 


where 


dk,mn  ^k+l,mn  qk,mn  ^k+l,mn  qk,mn  ~  nk,mn 
The  unnormalized  mean  squared  error  in  the  object  domain  is  defined  as 


E0k  ~  iC  ^k,mn^  iLrf^k.mn^ 
mney  nn 

And  by  applying  Parseval's  theorem,  Eq.  (A2)  becomes 

r2  /„  \2 


EFk  -E  (nk,mn>' 


In  the  following  we  will  show  that  |  mn  |  £  |  n|<  mn  I  f°r  (m,n), 

proving  that  Egk  <_  Epk.  (1)  For  (m,n)  where  g'k  mn  satisfies 
the  constraints  M.e.,  fm.n)  Y],  dk>mn  =  0,  and  so  \dk'm\  <  |nk^mnJ 

(2)  For  (m,n)  where  a’k  mp  is  negative,  gk+-|  mn  =  0  and,  from  Eq.  (A7) 


dk  ,mn  ^k  ,mn 

”^k,mn  “  nk,mn 

Since 

^  k,mn 

is 

neoative,  \ 

is  nonnegative 

^k  ,mn 

is  nonneqative  we  have 

dk,mn  -  dk,mn  * 

%,mn  =  ~  nk,mn 

and  sc 

1  |  dk ,mn 

1 1 

1  nk ,mn  I-  0)  For 

(m,n)  where  the 

zero  (for  the  case  in  which  the  object's  diameter  is  known  and  that  con¬ 
straint  is  imposed),  gk+1>mn  =  gk>rnn  =  0,  also  leaving  d^  =  - 
nk ,mn’  and  so  ldk ,mn  I  <l"k,mnl  *  Therefore,  since  |dk>fJ  < 


n.  I  for  all  (m,nl,  then  from  Eqs.  (A8)  and  (A9)  we  have 
iv  9  mn  I 


E^  < 

t0k  ^  tFk 


Fourier  transforming  Eq.  (A7), 


\,pq  *  Gk+1  ,pq  Gk,pq 


(A10) 


(All) 
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and  so 


!°k,pc,!2  '  l<Vl,pql?  *  l^.pql2  -  2  ReCG^pq  Gk*l,pq]  "  <A,2> 

where  an  overbar  denotes  complex  conjugation  and  Re  the  real  part.  From  Eq. 
(A3), 

Nk+l,pq  =  Gk+1 ,pq  "  Gk+l,pq 

and  so 

iVl.pql2  -  lGk+l  ,pql  ^  *  IVl.pql2  *  2  ^l.pq  Gk.l,pq]  <A14> 
Taking  the  difference  between  Eq.  (A12)  and  Eq.  (A14)  and  using  the  fact  that 
|G'k,pq  I  -|GVl,pql  “  I  Fpq  I  *  we  have 

l°k,pql  "  lNk+l,pq!  “  2  Re^Gk+1,pq  ^k+l.pq-^  “  2  Re^Gk,pq  ^k+l.pq-^ 


*  2lFpql*lGk+l,pql  C1  -  cos  ®pq]  >  0 
where  e(u)  .  argfG • k ^ pq  6^  pq].  And  therefore  for  all  (p,q), 

iNk+l,pq  I  ±  lDk,pq  I 


( A15 ) 


(A16) 


This  result  can  also  be  easily  seen  from  the  geometry  of  Figure  Al.  The  re¬ 
sult  Is  the  same  irrespective  of  whether  |  |  is  greater  than  or  less 

than  lFpql.  Therefore,  by  using  Parseval's  theorem  on  Eq.  (A8)  and  usina 
Eq.  (A2),  we  have 


EF,k+1  ^  E0k 

Combining  this  with  Eq.  (A10)  gives  the  desired  result 

2  2  2 
EF,k+l  i  E0k  ^  EFk 


(A17) 


( Al  8 ) 


That  Is,  the  unnormalized  mean-squared  error  must  decrease,  or  at  least  not 
Increase,  with  each  Iteration  of  the  error-reduction  approach.  The  quanti¬ 
ties  usually  considered  are  the  normalized  mean-squared  errors,  Eqs.  (4)  and 
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(5).  Since  these  are  simply  proportional  to  the  unnormalized  mean-squared 
errors,  Eq.  (A18)  also  holds  for  the  normalized  mean-squared  errors. 

Complete  staqnation  of  the  error  reduction  approach  could  occur  if 
dk  mn  =  -  nk  mn  for  all  (m,n),  in  which  case  from  Eq.  fA 7)  one  gets 
^k+l  mn  =  mn  ^or  Another  condition  of  stagnation  is  for 

e  =  0  for  all  (p,q)  in  Eq.  CA15). 
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